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ABSTRACT 

It has been previously shown that moons of extrasolar pla nets may be detectabl e 
with the Kepler Mission, for moon masses above ~ 0.2 M e (Kipping et al.l [20090). 
Transit timing effects have been formerly identified as a potent tool to this end, ex- 
ploiting the dynamics of the system. In this work, we explore the simulation of transit 
light curves of a planet plus a single moon including not only the transit timing effects 
but also the light curve signal of the moon itself. We introduce our new algorithm, 
LUNA, which produces transit light curves for both bodies, analytically accounting for 
shadow overlaps, stellar limb darkening and planet-moon dynamical motion. By build- 
ing the dynamics into the core of LUNA, the routine automatically accounts for transit 
timing/duration variations and ingress/egress asymmetries for not only the planet, 
but also the moon. 

We then generate some artificial data for two feasibly detectable hypothetical systems 
of interest: a i) prograde and ii) retrograde Earth- like moon around a habitable-zone 
Neptune for a M-dwarf system. We fit the hypothetical systems using LUNA and demon- 
strate the feasibility of detecting these cases with Kepler photometry. 

Key words: techniques: photometric — planets and satellites: general — planetary 
systems — eclipses — methods: analytical 



1 INTRODUCTION 

In recent years, the possibility of detecting the moons of 
extrasolar planets , so-called "exomoons", has received in- 
creased attentio n (Sartoretti & Schnci derll 19991 ; lHan fc Hani 
20021: ISzabo et al.l|2006j;ISimon et al.ll2009l i 



- .Lewis et al.ll2008l : 

Kipping! l2009al lbl; lsato fc Asadall2009l : lKippinell2010bl ). Ex- 



trasolar moons may be frequent, temperate abodes for life 
and a determination of their prevalence would mould our 
understanding of the abundance of life in the Universe. Al- 
though many techniques have been proposed, it is the tran- 
sit method which seems to offer the greatest potential to de- 
tect h abitable-zone moons in the near- future |Kipping et al.l 
l2009d ). 

In general, there are two ways in which a moon can 
be identified using transi ts. The first of these is the m oon's 
own transit light curve (|Sartoretti fc Schneider|[l999l ) and 
the second is the family of techniques known as "transit 
timing e ffects" . This includes both tr ansit timing variations 
(TTV) jSartoretti fc Schneiderlll999l ) and transit duration 
variations (TDV). Further, TDV has two different compo- 
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nents- one due to velocity variations (TDV-V) (|Kippinel 
2009aH and one due to trans it impact parameter variations 
(TDV- TIP) (|Kippindl2009bl ). TTV and TDV-TIP are both 
due to the position of the planet varying in response to the 
moon's presence (positional wobble). In contrast, TDV-V is 
due to velocity variation of the host planet in response to the 
moon's presence (velocity wobble). Co mbining all of th ese 
effects allows for a unique solution (see [Kipping 2009bl for 
details) for the moon-to-planet mass ratio and exomoon pe- 
riod (and therefore orbital semi-major axis through Kepler's 
Third Law). 

To perform the observations, highly precise, contin- 
uous photometry is required and the Kepler Mission is 
thought to be the best instrument currently available to 
conduct such a search (|Kipping et al.l l2~009d ). So far, two 
searches for exomoons have been conducted using the Ke- 
pler photometry, utilizi ng the transit timing te chniques on 
Kep ler-4b through 8b (IKipping fc Bakosll2011al ) and TrES- 
2b (|Kipping fc Bakosl l2011bl ). Despite null- results (which 
was expected as all of these planets are hot-Jupiters), the 
studies indicate sensitivity into the su b-Earth mass regime, 
as predicted by IKipping et al.l (|2009d ) . The transit timing 
equations are all built around the premise of a constant 
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planet-star separation and constant planetary velocity dur- 
ing the timescale of the transit. Whilst this is generally a 
good approximation except for moons on very short peri- 
ods, there is an obvious desideratum to make the expressions 
dynamic. 

A possible problem with the current technique is that if 
the moon's light curve is detectable directly, then not only 
is this an extra piece of data we are ignoring but also it 
could fundamentally invalidate many of the transit timing 
methods. The ideal solution would be to simulate both the 
planet and moon transit light curves, including limb darken- 
ing which is critical for Kepler. In addition, these light curves 
could be computed with the TTV and TDV effects inher- 
ently built into the model in a dynamic way. Not only would 
this improve upon the static approximations made in the 
TTV/TDV equations, but it would also bring in TTV/TDV 
of the moon itself. 

Before constructing such an algorithm, one must keep 
the purpose of such a routine in mind. This algorithm will 
be used to search for exomoons which will be done by fitting 
light curves with the new algorithm. Given the large number 
of free parameters and inevitable intricate inter-parameter 
dependencies, a Monte Carlo based method seems required. 
Since such methods are inherently computationally expen- 
sive making often millions of calls to the simulation rou- 
tine, then it is clear that any algorithm we design must be 
very fast to execute. This essentially excludes methods based 
upon pixelating the star or other numerical methods. The 
clear requirement then, is a completely analytic algorithm. 
The list of requirements are: 

■ Analytic (absolutely no numerical components) 

■ Dynamic (inherently accounts for all timing effects) 

■ Limb darkening incorporated (including non-linear 
laws) 

■ All orbital elements accounted for (e.g. eccentricity, lon- 
gitude of the ascending node, etc) 

In summary, such an approach would not only offer 
many advantages over the previous timing methods, but 
would also be highly practical for conducting a batch-style 
search for exomoons in archival data. In this paper, we in- 
troduce the fundamental equations needed to construct this 
algorithm, known as LUNA, and apply it in some hypothetical 
examples. 

In ij2j we describe the framework for computing the 
transit light curves of a planet with a moon and outline the 
assumptions made. Accompanying details on the model used 
for the sky-projected motions of the planet and moon can be 
found in the Appendix ( §A|) . In Ej3j we provide expressions 
for the eclipsing area of the moon in front of the star in all 
27 possible principal configurations. In fj4] example transits 
light curves are generated and re-fitted for habitable-zone 
exomoons detectable with _ftTep/er-class photometry. Finally, 
we discuss comparison to previous methods proposed in the 
literature and provide concluding remarks in [JS] 



2 LIGHT CURVE GENERATION 

2.1 Converting Time to True Anomaly 

To generate a light curve, the relative positions of the star, 
planet and moon must be calculated at every time stamp in 



the photometric data set in question. In particular, one re- 
quires the sky-projected separation between the planet and 
star, Sp* , the moon and the star, Ss* and the planet and the 
moon, Sps- We direct the reader to the Appendix for 
details on the derivation of the expressions for these terms, 
where we employ an analytic approximation for the three- 
body problem w hich we dub as the nested two-body model 
l|Kippindl20i0bh . The coordinate system employed through- 
out is also defined in the Appendix fS|A")l. It can be shown 
that these three S terms are a function of fs* and fsB which 
may be written, in general, as a function of time alone. Note 
that we define Jb* as the true anomaly of the planet-moon 
barycentre around the star and fsB as the true anomaly 
of the satellite around the planet-moon barycentre. There- 
fore, a necessary prerequisite is to convert time stamps to 
true anomalies. Let us begin with considering the standard 
practice for a planet by itself. 

2.1.1 The Planet- Only Case 

The planet-only case is a very familiar one. However, we 
will cover the conversion of time to true anomaly carefully. 
Writing out the steps explicity will allow us to identify the 
necessary procedure for the more complicated planet-moon 
case we will soon face. 

Typically, we have a time series running from an initial 
time stamp of ti up to tjv (where N is the total number 
of measurements) with the time of transit minimum occur- 
ring at r. Our first task is to convert the time array into 
an / array. To accomplish this we must define a reference 
point in time at which the true anomaly is known. Whilst 
we are free to use any reference point we so desire, a typ- 
ical choice is the time of transit minimum, r. In the ex- 
ample above, this choice it will change our array to run 
from t'i(= ti — t) — > t' N (= tjv — t) with the transit cen- 
tred at time t' = 0. The true anomaly at the time of transit 
minimum can be found by solving dS/df = for /. For 
even a planet-onl y case, this is non-trivial and leads to a bi- 
quartic equation l|Kippin g]|2008l ). The solution may be found 
by a series expansion about the point of inferior conjunction 
fit — t) — 7r/2 — uj + n, where rj repre sents a perturb ation 
term, which is given up to 6 th -order in lKippind (|201ll ). 

The third and final step is the application of Kepler's 
Equation. With the reference true anomaly known, this is 
conve rted to a reference mean anomaly via the usual rela- 
tions (|Murrav fc Dermott|[l999h . Then, the mean anomalies 
at all times may be calculated since this parameter scales lin- 
early with time. Finally, the mean anomalies are converted 
into true anomalies using Kepler's Equation. 

So to summarize we have three steps: i) subtract the 
time of transit minimum from all times ii) define the ref- 
erence anomaly as the time of transit minimum iii) assign 
mean anomalies for all times, which are then converted to 
true anomalies using Kepler's Equation. 

For data spanning multiple transits, the reference time 
should be close to the weighted mean transit number. This 
selection typically minimizes the correlation between the or- 
bital period and reference time. In such a case, the range of 
allowed values for the fitting routine to explore vary from 
t = (r gucss — P/2) — > (r gucss + P/2). Moving outside of this 
range would cause the fitting routine to assign a different 
transit epoch as the reference transit instead. 
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2.1.2 The Planet-Moon Case 

Now consider a planet with a moon. The barycentre of the 
planet-moon system essentially behaves in the same way 
which was used to describe the planet-only case. The ob- 
served time of transit minimum of the planet is displaced 
from a linear eph emeris by a sm all time St, due to transit 
timing variations (|Kippinj2009ah . However, the barycentre 
behaves identically to before and still passes the star at a 
minimum when fs* = (w/2 — usb* + ff). Therefore, there is 
no need to change anything here except that it is understood 
that the r value we fit for is not the time of transit mini- 
mum of the planet but the time of transit minimum of the 
planet-moon barycentre across the star, and thus we denote 
it as tb*. 

For the moon, an analogous logic may be followed. 
There must exist a second time shift, tsb, to account for the 
phasing between the moon in its orbit around the barycen- 
tre. It also seems clear that this phase time has a range of 
Psb (orbital period of the moon around the barycentre), in 
analogy to the tb* case. Let us imagine we have subtracted 
tb* from all of our times so are left with times running from, 
say, t[ — > t' N . The obvious deduction is that we must now 
subtract a value tsb to get to a new time frame t", just as 
we did for tb*. 

A natural choice for tsb is the instant when 
dSsB /d/ss = 0. In analogy to the previous case, this occurs 
at fsB — tt/2 — lo sb + rjsB, where t)sb is the perturbation 
term. Because this parameter is a phasing term, we prefer 
to use 4>sb = {2tvtsb)/Psb as a definition, to decrease the 
correlations with the orbital period. 

2.2 Retrograde Orbits 

A unique problem with LUNA is that exomoons can be either 
prograde or retrograde. Whilst for planets this is also true, 
it actually makes no difference to the transit light curve (al- 
though the Rossiter-McLaughlin phenomenon is affected). 
However, for a moon the sense of motion is di stinguishable 
via the TTV and TDV effects (|Kippindl20TObT ), meaning it 
does affect the transits and so cannot be neglected. 

We tackle this by treating retrograde moons as having 
7r < isB < 2n and prograde m oons as having < isB < 7T. 
It is shown in iKippind (|201ll ) that this definition produces 
the correct retrograde behaviour for the selected coordinate 
systerrQ. 

We point out that the asymmetry between prograde 
and retrograde moons is small. The source of the asymmetry 
comes from the relative phase difference between the TDV- 
V and TDV- TIP effects, which i s for prograde moons and 
7r for retrograde (|Kippindl2009br i . Therefore, the determina- 
tion of the sense of orbital motion is only generally possible 
if both TDV-V and TDV-TIP effects are observable. 

2.3 Small-Moon Approximation 

From a naive perspective, one might start by considering 
that the transit light curve of a planet and a moon can be 

1 Note that the longitude of the ascending node of the moon 
in its orbit around the planet-moon barycentre has the range 
s? fl S B < 2tt 



computed by calculating the light curve of both separately 
and then simply adding the signals together. However, there 
exists numerous scenarios where this would break down. For 
example, the moon could be eclipsing the star but fully in- 
side the planetary disc and thus the change in flux caused 
by the moon would be zero. 

To overcome this, we start by generating the planetary 
light curve in the nor mal way. This can be done using the 
IMandel fc Ago! (|2002l ) routine with any limb darkening law 
we wish and not making any size approximations. Having 
computed this curve, we need to add on the contribution 
from the moon. What really matters is what part of the 
moon is "actively" transiting the star. We define this as 
the area of the moon which overlaps the star but does not 
overlap the planet. If we can find this area, then we could 
compute a light curve for a planet + moon without any limb 
darkening effects immediately, which is a necessary first step. 
The subject of the moon's actively transiting area will be 
derived later in Sj3] 

For now, let us assume we know what this area is and 
continue to think about how the limb darkened light curve of 
the actively transiting moon can be computed. In general, we 
are interested in achieving two things i) focussing on moons 
rather than, say, binary-planetfl ii) actually detecting such 
signals. The first of these means we are dealing with small 
objects such that s < 0.1 in all cases and most likely s < 0.01 
in most cases (where s is the ratio of the satellite-to-star 
radii). The second point means that we need to be able 
to perform fits of transit light curves including all of the 
planet and moon properties. The inevitably large amount of 
parameter space necessitates computationally efficient and 
expedient algorithms. Putting these two arguments together 
indicates that the best way forward is to employ th e smal l- 
planet approximation case used in IMandel fc Agoll (|2002h . 
This will ensure accurate modelling of the limb darkening 
but very fast algorithms. 

We stress that the planet is still modelled using the 
full equations since the host planet could be a Jupiter sized 
object, yielding p > 0.1 (where p is the ratio of the planet-to- 
star radii). Therefore, we assume the moon is small and the 
planet is not. One limitation of this assumption is "binary- 
Jupiters" . If we have a binary-planet system for which the 
smaller body satisfies s > 0.1 then the algorithm we develop 
here will be limited in accuracy. However, the current focus 
of LUNA is solely to detect small bodies. 

2.4 Limb Darkened Light Curve for the Actively 
Transiting Lunar Component 

2.4-1 Uniform source 

Let us assume that the component of the moon's area which 
actively transits the star is known and equals As, transit- Be- 
fore we can compute the resulting limb darkened light curve 
from this component, we must first consider the case for a 
uniform source. 

We start by assuming the actively transiting component 
of the moon is equal to that which would transit in the ab- 
sence of a planet. In other words, we here ignore the planet 

2 Although we intend to extend LUNA to binary-planets in future- 
work. 
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and will re-introduce i t later. We follow the methodology of 
iMandel fc Ago] i|2002h . where the ratio of obscured to unob- 
scured flux is given as F§(s, Ss*) = 1 — Ss*) (replacing 
the appropriate symbols for the exomoon case): 





i [s 2 K0 + Kl — K2] 

s 2 

1 



1 + s < S s , 
\l-s\ < Ss* 
S a ,^l-a 
S s , < s- 1, 



1+8 



KQ = COS 



(1) 
2 + 



where «i = cos^ 1 [(l - s 2 + Sf, )/2Ss 

Sf. - l)/2sSs*\ and k 2 = 0.5^/451, - (1 + 5|, - s 2 ) 2 . 

The most interesting case is clearly (1 — s) sC Ss* ^ 
(1 + s), where only a portion of the moon transits the star. 
The area of the moon which transits the star may be denoted 
as as* where the subscript indicates S is transiting *. 

In this case, A| = as* /it, but a really generally de- 
scribes the area of intersection between any two circles. We 
prefer to generalize the equations at this stage, which will 
prove useful later. Therefore, the area of intersection caused 
by object of radius r transiting object of radius R, with 
separation S, is: 



Table 1. List of cases identified by IMandel fc Agoll ^2002h . We 

use the same case notation in this work, but modifying them to 
emphasise the focus on moons. 



Case 


Condition 


Area of star which is eclipsed 


I 


1 + s < S s * < 00 





II 


l-s<S s *<l + s 


as* 


III 


s < S s , < 1 - s 


ns 2 


IX 


< S s , < s 


■KS 2 



a(R,r, S) = r 2 ti (R,r, S)+R 2 K 1 (R,r, S) - n 2 (R,r, S) 

(2) 

/ p on \ S 2 + r 2 -R 2 ] 
Ko{R, r, S) = arccos ^ — J (3) 

I S 2 + R 2 — r 2 1 
ki (R, r, S) = arccos [ J (4) 

m / 4S 2 R 2 -(R 2 + S 2 -r 2 ) 2 
K 2 {R,r,S) = ^ (5) 

Let us now re-introduce the planet. The consequences 
are that in various circumstances the portion of the moon's 
shadow which actively transits the star is diminished due to 
overlap between the planet and the moon. We will discuss 
the derivation of the actively transiting lunar area in §j3j but 
for now it is sufficient to say that for a uniform source the 
loss of light due to the lunar component is: 



rize the different cases of interest, we direct the reader to 
Tabled] 



Fl(s,S s ,) = l- As ' tlansit (6) 

As an example, in the case of no overlap between the 
planet and the moon, but the moon fully inside the stellar 
disc, As jt ransit = its 2 and thus we recover the expected form 
shown in Equation [T] 



2.4.2 List of Cases 

Looking at the equations for a uniform source (Equation [lj, 
one may identify three distinct cases: i) moon fully outside 
the stellar disc ii) moon partially inside the stellar disc iii) 
moon fully inside the stellar disc. The first of these is triv- 
ial; there is no transit occurring. The third is also trivial for 
a uniform source but not so for a limb darkened one. For 
reasons which become clearer later, it is necessary to sepa- 
rate ii) into two different cases dependent upon whether the 
moon overl aps the stellar centre o r not (cases III and IX in 
the orig inal IMandel fc Agoll ( 2002) assignment). To summa- 



2.4.3 Case III 

Under the assumption of a small ratio-of-radii (s <^ 1) one 
may assume the surface brightness of the star is constan t 
under the disc of the transiting body l|Mandel fc A gol 2002). 
We will here assume non-linear limb darkening of the form 
I(r) = 1 — S* =1 Cn(l — M n//2 )> where I(r) represents the 
radial intensity distribution from the star with 1(0) = 1, fj, = 
cos 8 = y/1 — r 2 and r is the normalized radial coordinate 
on the disc o f the star. The non-linear limb darkening law 
l|Claretll2000l l is chosen as it may be easily used to compute 
lower order limb darkening laws, specifically quadratic and 
linear laws. 

Let us draw an annulus inside the star with width 2s at 
radius Ss* to represent this surface. We begin by considering 
case III, where the body is fully inside the star but does not 
cover the stellar centre (i.e. s < Ss* < 1 — s). We have: 
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f 2rl(r) 
Jo 



dr 



^ n + 4 



F. 



S.annulus 



(Ss*,s) 



n=l 
S S «-s 



rCi 



(7) 



2r/(r) dr 



( 1 — bm — V^l — a-mj 
+ (5|,-1) (VT - 6^- #l-a m ) 



2 

+ 21 



+ sSs* (^2^1-a m +2v / l-6 m -5 

7C 2 |^(S 2 ^Vl - 6m - VI - Urn) 

+ (Si, - 1) (V 1 - fem - VI - a m ) 

+ 2sS s » (Vl-am + Vl-bm-3)] 

-6[c 3 ( S 2 ((l-am) 3/4 -(l-6m) 3/4 ) 

+ Sl,((l-a m ) 3/4 -(l-6m) 3/4 ) 

+ s5 s , (-2 (1 - a m ) 3/4 - 2(1 - & m ) 3/4 + 7) 



(l-a m ) 3/4 + (l-bm) 3/4 ) 

0] 



+ 7sS s * (c4 (s 2 + S 2 S , 



(8) 



where a m — (Ss* — s) 2 and b m = (Ss* + s) 2 and the m 
subscr ipt denotes this definition comes from lMandel fc Agoll 
(2002). A further simplification is possible by using a mr = 
(l-a m ) 1/4 and b mr = (1 - 6 m ) 1/4 . 



aIU transit 

A S - -7m 

S, annul us 



AT 



tv[(Ss, + s) 2 - (S s , - s)2] 

-^S.transit 



7r(fo m - «m) 

The final flux is given by: 



FS = 1- M 



S'.an mil us 
-ftotal 



(10) 



(11) 



(12) 



244 Case IX 

If the transiting body has an equatorial orbit, such that 
Ss* < s, the setup for case III is invalid since one of the 
integration limits is negative and r is strictly defined to be 
r ^ 0. This gives us case IX: 



1 S, annulus 



S S „ + s 



1 

+ 14 



2r/(r) dr 



/ 5 

1 - femr + (S + S S *) 2 y ~ ~ + bmr 



1 - b 2 mr + (s + S s *) 2 ( - | + b mr 



U(s + S s *) 2 ~7{s + Ss*) 4 c 4 



Cl 



+ 2(4 - 4&L + (s + 5s,) 2 (-7 + 4&L-))c 3 

(13) 

For the ratio of the annulus to actively-transiting-lunar 
area we have: 



5,111 

S, annul 



4 f 2 2 
us {Ss*,s) = g [ s (b mr — a mr ) + (S s * — l)(b m r — Omr) 

+ 2sSs* ( ~ ^ + amr + k""-)] c i 



2 

+ 21 



IX -As, transit 



(14) 



7^S 2 (&mr - «mr) + (Ss* ~ 1) (&mr - « 
+ 2s5s* (-3 + a 2 „ r + &mr)) c 2 

^ y ^mr a mr + 2s Ss * (~ a mr & mr ) 

+ s 2 (a mr — b' mr ) + Ss*(a mr — b mr )jC3 



+ 7sSs*{-l + (s 2 + S 2 s,)a) 



(9) 



The final step is to acknowledge that the transiting 
body does not transit an area equal to that of the annulus, 
but in general, a smaller area. Therefore, we must multiply 
the flux from the annulus by a ratio-of-areas given by: 



24.5 Case II 

The final case we need is case II. This is when the body is 
in the ingress/egress portion i.e. (1 — s) < Ss* < (1 + s). In 
this case the annulus flux is given by: 



S, annulus 



2rl(r) dr 



Ss*- S 



-(-5 + 4a mr ) ^tnr '-'1 



5 

a m — 1 
42 



- 42 + (42 - 28a mr )e 2 + 6(7 - 40^)03 



+ 21(1 + a m )c 4 



(15) 



The ratio of the annulus to actively-transiting-lunar 
area is: 
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^11 &S t transit 

7r(l — a m ) 



2.5 Summary 

The light curve of the actively-transiting-lunar compo- 
nent has now been evaluated, but we have not yet de- 
rived As itrans i t : the sky-projected area of this component. 
The planetary l ight curve is computed using the usual 
IMandel fc Agoll (|2002h routine, utilizing Sp*. The only re- 
maining task is therefore to find As itra nsit, which we deal 
with next in i|3] 



3 ACTIVELY TRANSITING LUNAR AREA 
3.1 Principal Cases 

The one-body case is described by one parameter, S, which 
can have three states. As we saw in ^J2j the light curve of 
the actively-transiting-lunar component can be generated in 
each case by just knowing Ss* and A^transit, in the small 
moon approximation. In what follows, we always assume 
s 2 < p 2 < 1 where s = Rs/R*. It is therefore clear that 
our task is to find A^transit in all possible configurations i.e. 
the actively transiting area. This will enable us to compute 
A e in cases II, III and IX and thus produce limb darkened 
transit light curves for the moon. 

With three S values now in play (Sp*, Ss* & Ssp) 
more cases are possible than before. In total 3 3 = 27 prin- 
cipal cases exist, but many of these are unphysical under 
the conditions that R > p > s (where R = 1) and that 
Sps < \Sp* — Ss*\. In Table [21 we list all the possible cases 
by permuting each S parameter into each of their respective 
three states. The cases which are unphysical are marked ac- 
cordingly. Figures [T] [2] & illustrate example configurations 
for each principal case. 

Table [5] shows that for some cases, the amount of lunar 
area actively transiting the star is dependent upon the posi- 
tion of the planet (these case are marked with a *). This is 
because the planet's shadow overlaps with the lunar shadow 
and so the the moon cannot block out to its full potential. 
We label these cases as exhibiting "areal interaction". 



3.1.1 Case 26 

We start by considering the simplest of the areal interac- 
tion cases, case 26 (see Figure [3]). We denote this as £26 in 
mathematical notation (and likewise for the other case num- 
bers). £ is chosen in a nalogy to the case notation e used by 
IMandel fc Agoll | |2002j) . For case 26, the planet and moon 
are both completely inside the stellar disc but are inter- 
acting with each other. This is very similar to the case II 
situation for a one-body transit where the star has become 
the planet and the transiting-body has become the moon. 
So {R, p, S} — > {p, s, Sps}- The area of the moon which left 
actively transiting the star is given by ^4f 2 t^ anB it = 7rs 2 — asp. 



3.1.2 Case 17 

The next most complicated case is case 17. Here the planet 
is on the limb of the star and the moon is on the limb of the 
planet, but completely within the star. It can be seen that 
^transit = -4f 2 t 6 ran S it = tts 2 - a S p and thus we already have 
the required equation in hand. 

3.1.3 Case 23 

The third most complicated case is when the planet is com- 
pletely inside the star but the moon sits on the planetary 
limb and coincident ally the stellar limb. If the planet was 
not present, the moon would eclipse an area given by as*. 

However, the planet is present, and it blocks out a por- 
tion of this area. Under the conditions of case 23, the over- 
lapping area between the planet and moon must always lie 
completely within the stellar disc since the planet always 
lies completely within the stellar disc. It can be therefore be 
seen that the overlapping area, asp, will always be less than 
as*. Therefore, the final area which is actively transiting the 
star from the moon alone is: 

^"transit = «S. ~ «SP (17) 

3.2 The Sub-Cases of Case 14 

Case 14 is the most complicated case to consider. Indeed, 
case 14 actually has multiple sub-cases which define the var- 
ious possible behaviours. Although similar to case 23, the 
planet's shadow now does not completely eclipse the star 
and so the overlapping area between the moon and planet 
also does not necessarily have to lie completely within the 
stellar disc. As a consequence, multiple different states exist 
and a singular expression is not possible for -As.translt- 

The first key question is: how many sub-cases actually 
exist? The problem we are thinking of is that of three cir- 
cles interacting to give an overlapping region, for which we 
require the associated area. The problem of finding the area 
of overlap between th ree circ l es wa s studied extensively in 
the pioneering work o f iFewelll l| 2006). Indeed, an analytic so- 
lution for the overlap of three circles did no t exist until 2006 
when Fewell presented the solution. IFewelll (|2006T ) identified 
nine possible cases in total, which w e label as Few ell case 
(1) -> (9) and viewable in Figure 5 of lFewelll l|2006h . To as- 
sign identities in our framework, we label the biggest circle 
the star, the middle-sized circl e the p l anet a nd the smallest 
circle the moon in Figure 5 of Fewell (2006). 

Going through the IFewelll (|2006fl scenarios with our as- 
signed idejitities 1 _the first thing we notice is that only four 
of the IFewelll | 2006) cases remain consistent with the case 
14 conditions, the others are much simpler scenarios which 
we dealt with earlier as the principal cases. I n Tabl e El we 
show the corresponding case numbers to the IFewelll (20061 ) 
cases for completion. 

The|Fewcll (2006]) solution for the area of common over- 
lap between three circles is based upon a formula for the 
circular triangle, which uses the chord lengths and radii of 
the three circles involved. Only in Fewell case (1) doe s a cir- 
cular triangle act ually e x ist an d this is the focus of I Fewell 
l|2006h . However, IFewelll | 2006) performs the derivation by 
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Table 2. List of principal cases. For each case we label whether the case is geometrically physical or not. We also give the expected 
eclipsed area due to each body. More complicated cases are marked with a * and exhibit arcal interaction. ** indicates multiple sub-cases 
exist, which is discussed in 13.21 



CbjSC Number [£*] 




S P 














Sps 


PliysicS/l? 


^P,transit 


4di 

-^■0, transit 


i 




S P , > 1 


+ p 






S S * 3= 1 + s 






Sps 3= P + s 










2 




Sp, ^ 1 


+p 






Ss, > 1 + s 




p — 


s < Sps < P + s 










3 




Sp, ^ 1 


+p 






S s « 3= 1 + s 






Sps ^ P — s 










4 




S P , ^ 1 


+ p 




1 


- s < S s , < 1 - 


h s 




Sps 3= P + s 







"S* 


5 




Sp, > 1 


+ p 




1 


- s < S s » < 1 - 


hs 


p- 


s < Sps < P + s 







"S* 


6 




Sp, ^ 1 


+p 




1 


— s < 05* < 1 - 


h S 




Sps =S p - s 


X 






7 




Sp, ^ 1 


+p 












Sps 3= P + s 


/ 
V 







8 




Sp, ^ 1 


+p 






OS* ^ 1 — S 




p - 


s<Sps<p + s 


X 






9 




Sp, ^ 1 


+p 






Ss* ^ 1 — s 






Sps ^ P - s 


X 






10 


1 


- p < Sp, < 1 - 


\-p 




S s * 5= 1 + s 






Sps 3= P + s 




Q:p* 





11 


1 


-p < Sp, < 1 - 


hp 




S s * 3= 1 + s 




p — 


s < Sps < p + s 




a P* 





12 


1 


-p< Sp, < 1 - 


\-p 




S s , 3= 1 + s 






Sps ^ P — s 




C£ p* 





13 


1 


-p < Sp, < 1 - 


hp 


1 


- s < S s , < 1 - 


Is 




Sps 5 s P + s 




ap* 


as* 


14** 


1 


-P < Sp, < 1 - 


hp 


1 


- s < S s , < 1 - 


h s 


p- 


s < Sps < P + s 


v 7 


&P* 


variable 


15 


1 


-P < Sp, < 1 - 


hp 


1 


- s < S s . < 1 - 


h s 




Sps <P~ s 


V 


&P* 





16 


1 


- p < Sp, < 1 - 


hp 




Ss, s: l - s 






Sps 3= P + s 


V 


a P* 


7TS 2 


17* 


1 


-p < Sp, < 1 - 


hp 




S s * < 1 - s 




p- 


s < Sps < p + s 


V 


a P* 


7TS 2 — CUSP 


18 


1 


-p < Sp, < 1 - 


hp 




Ss. < 1 - s 






Sps < P - s 


V 


ctp* 





19 




Sp, s= 1 


-p 






S s , 3= 1 + s 






Sps 3= P + s 


V 








20 




Sp, s= 1 


-p 






S s ,^l + s 




p- 


s < Sps < p + s 


X 






21 




Sp* < i 


-p 






S s , > 1 + 8 






Sps < P - s 


X 






22 




Sp, s: 1 


-p 




1 


- s < S s , < 1 - 


h s 




Sps 3= p + s 


V 


2 


as* 


23* 




Sp, sc 1 


-p 




1 


- s < S s , < 1 - 


Is 


p- 


s < Sps < p + s 


V 


2 


as* — asp 


24 




Sp, sc 1 


- p 




1 


- s < S s , < 1 - 


h s 




Sps ^ P - s 


X 






25 




Sp, sc 1 


-p 






Ss. < 1 - » 






Sps 3= P + s 


V 


2 

717T 


7TS 2 


26* 




Sp„ s= 1 


-p 






Ss* < l - s 




p- 


s < Sps < p + s 


V 


2 

7T/J Z 


7TS 2 — agp 


27 




Sp, sc 1 


-p 






Ss. < 1 - s 






Sps < P - s 


V 


2 






Table 3. List of lFewelll l|2006l ) cases and the corresponding case 
definitions used in this work. 



Fcwcll (2006) Case Number Case Number (This work) [£] 



(1) 


14.1 


(2) 


14.2 


(3) 


14.3 


(4) 


17 


(5) 


18 


(6) 


27 


(7) 


14.7 


(8) 


11 


(9) 


10 



X12 

yi2 



Intersection 12 

1 - p 2 + S'p, 



2Sp, 



2Sp, 



2S P ,(1 + p 2 ) - (1 - p 2 ) 2 - S%, 



(18) 
(19) 



calculating analytic expressions for the vertices of all three Intersection 13 

circles which may then be converted to chord lengths and v v i Q i . i , nn *. 

J ^tl3 = <-t 1 3 COS (7 — sin (7 [ZU ) 

these equations are useful for distinguishing between all the , , , ' , 

various sub-cases. ^ = sin 9 + ^ cos 9 ( 21 ) 

Let us consider three circles labelled 1, 2 and 3 of radii ; 1 — s 2 + S|, /22") 

7ii and positions {A'i, } and intersection points {Aij, 13 25s, 

(where two intersections points always exist for each i-j pair- , —1 / ~ ~ " ~ ~ 

ing under the case 14 conditions). Let us define {A^i} = = 2SsTV 2S s*( 1 + s ) ~ I 1 ~ s ) ~ S S* ( 23 ) 

{0, 0} and assume 1Z\ > IZ2 > Therefore, we have circle ^ S'p + 5| — Sp s 

1 being the star, circle 2 being the planet and circle 3 being cos 8 = 2 S'p Ss 

the m oon. The intersection points are given by (see iFewell 
(2006) for details of the relevant derivation): 



sine' = VT - cos 2 6' (25) 
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3.2.1 Fewell Case (1) 



Intersection 23 

X 23 = X23 cos 0" - y'23 sin 6" + S P , (26) 
J23 = X23 sin 6" + y 2 ' 3 cos 8" (27) 



What makes Fewell case (1) unique mathematically speak- 



„ _ p 2 — s z + S PS ,~g\ ing? i-e- under what conditions are the equations for (1) 

23 ~ 2Sp1; { ' valid? iFewelll |2006l ) provides expre ssions f or ev aluating; 

whether case (1) is satisfied or not. iFewe ll (2006) identi- 



^ 23 ~~ 2S PS \/ 2S Ps( p2 + s ^ ( p2 s ^ 2 ( 29 ) fies "two" conditions but one of these conditions is just the 

c-2 1 c.2 q2 case 14 conditions and need not be considered again. The 

nil _ °P* T '-'PS ~ &s* tm\ 

cos o — 2Sp Sp second condition is actually composed of two causal factors 

;; , and so in total there exists two conditions which uniquely 

sin 9 = yl — cos 2 6" (31) define case 14.1. Let us define conditions A and B as: 
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CASE 10 


CASE 11 


CASE 12 














CASE 13 


CASE 14 


CASE 15 




MULTIPLE SUB -CASES 


x s 


CASE 16 


CASE 17 


CASE 18 


^~C~3 


,' \ 





Figure 2. Principal cases 10 to 18. The star is black solid, the planet is black dashed and the moon is gray solid. 



Condition A 

{X12 - s s , cos e'f + (y 12 - s s , sin e'f < s 2 



(32) 



half the circle is involved and 14.1b occurs when more than 
half the circle is involved. The two cases are distinguished 
by evaluating the following: 



Condition B 

(#12 - s s , cos e'f + {y 12 + s s , sin e'f < s 2 (33) 

In order to have case 14.1, we require that condition 
A is satisfied and condition B is anti-satisfied (we call this 
condition B). The area atransit may be found by evaluating 
the circular triangle, principally defined by the chord lengths 
and radii. In fact, two sub-sub-cases exist here for 14.1 which 
we label as 14.1a and 14.1b. 14.1a is the case where less than 



Condition C 

Ss, sine' > y 13 + y J A - y v 13 (s s , cose' - x 13 ) (34) 

•-123 — <t-13 



If the above equation is true (condition C) then we have 



case 14.1a and ^overlap ma y De found using: 
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CASE 19 




CASE 22 



CASE 20 



CASE 21 



UNPHYSICAL 



UNPHYSICAL 



CASE 23 



CASE 24 



UNPHYSICAL 




Overlap = i\/(ci+C2+C3)(C2+C3-Cl) 
X sj (Ci + C 3 - C 2 )(C1 + C 2 - C 3 ) 

+ £ arcsin ^ - C ± ^1Z\ cl) (35) 
fc=i fc 

cl = (A- lfc -A- jfc ) 2 + (^ fc -J jfe ) 2 (36) 
Otherwise (condition (7), for case 14.1b, we must use: 



^fverlip = | V(Cl + C2 + C 3 )(c 2 + C 3 - Cl) 

X sj (Cl + C 3 - C 2 )(C1 + C 2 - C 3 ) 

3 

+ £(tt 2 * arcsin^-) 

fc — 1 

- j\J 4n i - c? - j\J 4ii 2 - % + j\f*n - c§ 

(37) 

The actively transiting area of the moon is now given 
by: ^"tansit = as. - Overlap- In the above equations, it 
should be noted that Ck represents the chord lengths. 
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3.2.2 Fewell Case (2) 

iFewelll l|2006h does not discuss in any detail how to compute 
the area of overlap for Fewell case (2), but this is simpler 
than Fewell case (1). In Fewell case (1) we had a circular 
triangle, whereas case (2) presents a circular quadrilateral. 
Before we consider the various scenarios, let us define under 
what conditions Fewell case (2) is satisfied. 

For Fewell case (2), the constraint is that the two inter- 
section points of circles 1 and 2 both lie outside of circle 3 
(this narrows it down to Fewell case (2) or (7)). Therefore 
we must satisfy conditions A and B. 

Next, the intersection points of circles 1 and 3 must lie 
inside circle 2 (this distinguishes it from Fewell case (7)). 
Assuming case 14, condition A and condition B, it is not 
physically possible to have a solution where one of these 
intersection points lies inside and the other outside of circle 
2. It is, however, possible to have both outside, which gives 
rise to Fewell case (7). Because of this fact, we only need to 
use one of the intersection points to define condition D: 

Condition D 

(X 13 -Sp.) 2 +yf 3 <v (38) 

But, even this is not en ough for a u nique solution. In 
Fewell case (2) of Figure 5 of lFewelll |2006), one could imag- 
ine mirroring the moon along the line connecting the two 
moon-star intersection points. This would still satisfy all of 
the above conditions but now no part of the moon is a ctively 
transi ting. The nominal shown in Figure 5 of lFewelj 

(2006), is labelled as 14.2a (shown in Figure|4)|, whereas the 
alternative case is 14.2b. Therefore, we define condition E, 
which if satisfied gives case 14.2a and if anti-satisfied gives 
case 14.2b: 



Table 4. List of sub-cases and conditions. A bar indicates the 
condition is anti-satisfied. 



Sub-Case 


Conditions Satisfied 


^S, transit 


14.1a 
14.1b 
14.2a 
14.2b 
14.3a 
14.3b 
14.7a 
14.7b 


A; B; C 
A; B; C 
A; B; D; E 
A; B; D; E 
A; B; F 
A; B; F 
A; B; D; G 
A; B; D; G 


overlap 
"S* ^overlap 

ns 2 — asp 


as* — ap* 
irp 2 — ap* - asp + as* 
as* 
as* - asp 


3.2.4 Fewell Case (7) 



For Fewell case (7), we must satisfy conditions A and B, as 
did Fewell case (2). The final condition for Fewell case (7) 
is simply the opposite of condition D, i.e. condition D. 

However, this is insufficient to give a unique solution. 
The standard Fewell case (7) has no interaction between 
the planet and moon within the stellar disc but this is not 
a condition, merely an artifact of the figure's construction. 
It is possible that the intersection points of 1 and 2 lies 
outside circle 3, and the intersection points of 1 and 3 lies 
outside 2, but still two distinct cases exist. These two sub- 
sub-cases are that if the intersection of 2 and 3 lies within 
1, then areal interactions must be occurring. We label this 
case 14.7b. Case 14.7a is that the circles do not interact. 

This first case occurs when both intersection points of 
2 and 3 lie outside of circle 1. It is not possible to have one 
inside and one outside under the previous conditions so far 
met. Thus we have condition G, which if satisfied yields case 
14.7b and if anti-satisfied gives 14.7a. 



Condition E 

(S s * -s)< (S P , - p) (39) 

Case 14.2a has Af^tansit — ' KS ' 2 ~ asp, otherwise for 
14.2b we have Af^it'^O. 



3.2.3 Fewell Case (3) 

To meet Fewell case (3), both intersection points of 1 and 2 
must lie inside circle 3 i.e. condition A and B. Sub-case 14.3 
contains two possible sub-sub-cases; 14.3a and 14.3b. The 
second sub-sub-case was identified quite late in our analysis 
and noticed during testing and debugging. 

The two cases are differentiated by whether the planet's 
centre is inside or outside the stellar disc. If outside, then 
we have case 14.3a, which satisifes condition F. 

Condition F 

S P , > 1 (40) 

In this case, 14.3a, ^| transit = a s* - ap.. For case 
14.3b, the planet is inside and we have -Af^r'ansit — T 2 — 
asp — ctp* + as*. 



Condition G 

xi 3 + yli < i (41) 

Accordingly, Ag 1 ^'^^ — as*. Otherwise the planet- 
moon intersections occur within the star and thereby neces- 
sitating a degree of interaction (case 14.7b, satisfying con- 
dition G), for which Af^nsit = «s* - a S p. 

3.2.5 Summary of Sub- Cases 

This completes every possible scenario for three circles in- 
teracting. The sub-case conditions are summarized in Ta- 
ble 0] and illustrated example configurations are shown in 
Figure 2] We have calculated the area of overlap in every 
case in an analytic manner. This will allow for extremely 
expedient computation of the light curve from a planet plus 
moon. 

Due to the plethora of cases, sub-cases and even sub- 
sub-cases, a large amount of testing of the LUNA algorithm 
has been executed. This was done by trying various system 
configurations and ensuring no gaps existed in the case con- 
ditions or unusual light curve features were produced by the 
algorithm. It was during this stage which we identified that 
sub-case 14.3 actually had two sub-sub-cases in the form of 
14.3a and 14.3b. 
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CASE 14.1a 




CASE 14.2 & 




CASE 14.1 b 



CASE 14.2 a 





CASE 14.3 a 



CASE 14.3 Z> 





CASE 14.7 a 



CASE 14.7 b 





Figure 4. Sub-cases of case 14- The subscript numbers (1,2,3,7) originate from the system employed bu \Fewell \2006l) . The ex t ra cas e 
conditions (a,b) come from additional complexity within each sub-case yet consistent within the original case definitions of Feweh (2006). 
The star is black solid, the planet is black dashed and the moon is gray solid. Relative sizes selected purely for depiction purposes. 



3.3 Flow Decision Diagram 



3.3.1 Decision 1 



With all the conditions stated, it is useful to construct a flow 
diagram of the decision route which LUNA should take. This 
flow diagram should be designed to minimize the number of 
calculations LUNA has to perform. We have conditions A to 
G for the sub-cases and three additional conditions from the 
principal cases, so the maximum number of decisions in any 
chain should be 10. Ideally, most decision routes will involve 
fewer than 10 decisions. 



The best-place to start is with the principal cases. The com- 
putation of the transit light curve will always require a com- 
putation of Sp*, as it is assumed that the planet always 
transits. 

Since we always have to compute this term anyway, the 
economic way to proceed is to base a decision tree upon its 
behaviour. This essentially gives us a three way-split into 
three decision trees: cases 1 through 9 (called "p-out"), 10 
through 19 ("p-part") and 20 through 27 ("p-in"). 
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3.3.2 Decision 2 

Having computed Sp*, should we compute Ss* and Sps? 
In several cases, only one of these is required to calculate 
the light curve. Consider, the p-out tree first (note that the 
second decision does not have to be the same in each tree). 
Since the planet is out of the stellar disc, areal interaction 
is impossible and thus Sps need never be computed. There- 
fore, the logical second decision for p-out is simply to eval- 
uate Ss* and see which of the three possible circumstances 
it is in. 

For the p-in tree, there are two possible questions one 
could pose. Remember, we wish to pose the question which 
leads to the most efficient algorithm. First possible question: 
evaluate Ss* and see which state it is in. Second possible 
question: evaluate Ssp and see which state it is in. 

If we ask the second question, the in-case requires no 
further computations. However, the part and out cases both 
require that we evaluate Ss* to proceed. 

If we ask the first question, the out case requires no 
further computation. However, the in and part cases require 
we calculate Ssp- Therefore, the two approaches seem to 
yield equivalent levels of complexity. In general, if the moon 
is transiting, Ss* will always be needed whereas Ssp is only 
needed if areal interaction occurs. Therefore, we prefer to 
ask about Ss* first. 

Finally, for the p-part tree, case 14 is the mix which is 
by the far the most complicated case computationally. We 
wish to choose a tree based upon minimizing our exposure to 
case 14. Ss* again makes a natural choice by the groupings 
of the actively transiting area formulas. 

3.3.3 Decision 3 

Since decision 2 is always the same, decision 3 is simply to 
evaluate Ssp and check its condition. There is no choice in 
this. 



3.3.4 Decision I4.I 

If case 14 occurs, we have a new decision tree opening up due 
to presence of sub-cases. The first decision is decision 14.1. 
Unlike the principal tree, each sub-case does not necessarily 
involve a calculation of all the possible conditions. There- 
fore, some conditions crop up more frequently than others. 
If all the conditions took the same time to compute, then we 
would prefer to compute those conditions which occur most 
frequently early on in our decision tree. 

Conditions A and B clearly provide the backbone to 
the decision tree and thus must be evaluated early on. Their 
evaluation is ultimately unavoidable so should be done at 
the start. Case 14.3 is unusual in that it can be identified as 
satisfying condition B only since all sub-cases anti-satisfy B. 
Condition B is therefore the obvious starting point in our 
decision tree. If verified, we have case 14.3 and if not we 
move into decision 14.2. 

3.3.5 Decision 14.2 

Since conditions A and B are the backbone and B has been 
evaluated, the next step is to evaluate A. If satisfied, we 



have case 14.1. Anti-satisfaction means we have either 14.2 
or 14.7. 

3.3.6 Decision I4.S 

If we have 14.1, then decision 14.3 must be to evaluate con- 
dition C. If we have 14.2 or 14.7, then the discriminator is 
condition D. 

3.3.7 Decision 14.4 

Not needed for 14.1 cases, but for Fewell cases (2) and (7) 
we must compute conditions E and G respectively. 

3.3.8 Summary of decision tree 

Figure [S] shows the decision process in a flow chart to sum- 
marize the approach of LUNA. 



4 EXAMPLE TRANSIT LIGHT CURVES 

4.1 Dealing with Inclination 

Transit modellers will be familiar with the pot-holes in the 
road presented when trying to fit transit light curves with 
the physical parameters. To start with, inclination is almost 
always in the range ~85 to 90 degrees and fitting for is* 
yields large correlations and lethargic routines. In practice, 
it is better to define the impact parameter of the transit 
chord across the star, bp*. 

b B , = {r B */R*)cosi B * (42) 

A similar strategy would seem advisable for moons and 
we can define an analogous quantity: 

bsp = (rsB/Rp)cosi S B (43) 

It should be noted that bsB is only the impact param- 
eter in the reference frame of the moving barycentre and 
will not be the observed impact parameter in the sky frame. 
Therefore, bsB can be greater than unity and yet the moon 
still transits the star. 

Another subtlety is that isB = 7r/2 + 8 would give a 
different light curve than isB = t/2 — 8, but both would 
yield the same value for bsB (in contrast to the situa- 
tion for a planet by itself). To account for this asymme- 
try, negative impact parameters are perceived as yielding 
isB — it — cos - [RpbsB /tsb]- Retrograde satellites with 
inclinations in the range tt < isB < 2-7T are adjusted in a 
second step where if a logical flag for retrograde orbits is 
switched on, n is added onto the inclination. 

4.2 Fitting Parameter Set 

A typical planetary parameter set would be 
{t, p 2 , {a/R*), b}. Due the str ong correlation between 
(a/R*) and 6 (|Carter et al.|[2008h . it is preferable to switch 
one of these parameters for a less correlated term. A typical 
choice is T, the duration for the planet to move from its 
centre overlapping the stellar limb to exiting under the 
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PRINCIPAL CASES 



PS 



O 

(/) 
CD 



£14. la 




SUB CASES 



Condition A 




Condition D 




Figure 5. Flow chart showing the decision process of choosing which case we are in. The decision tree has been optimized to marginalize 
the more CPU intensive computations until absolutely required. The principal cases are not shown but are easily seen in Tabled Square 
boxes indicate an evaluation of whatever is inside the box. Black arrows indicate a true statement and the dot-dashed gray arrows 
indicate a false. 
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same condition. An inver se-mappi n g expr ession to go from 
f to (a/P*) is given in iKippmel ^OlOal ). Note that a f 
equation for the moon is not possible since bsB > 1 is 
permitted and thus T would be complex. 

We now have the four fitted parameters 
{t~b*,P 2 , Tb*, &b*} for the planet-moon barycentre. 
Additionally, we need the orbital period, Pb», which may 
be fitted for if multiple transits exist. Radial velocity or 
secondary eclipse information may be used to expand 
the parameter set to include eb* cosljb* and e.B* sinwB*. 
Finally, the out-of-transit flux, OOT, and the blending 
factor, B (which cannot be fitted for) should also be 
included. The quadratic limb darkening coefficients ui and 
U2 give an additional two parameters bringing the total 
to 11. RV data may also require additional terms such as 
semi-amplitude, K, velocity offset, 7, etc. 

To include the exomoon, 9 new parameters are needed: 

{4>SB, (Ms/Mp), S, (asB/R*), b S B, PSB, es'BCOSOJSB, 

esB sinoiss, Q.sb}- In total, this means we have 20+RV pa- 
rameters. In practice, not all of these will be fitted for e.g. 
the blending factor. 

4.3 Example Simulations & Fits 

In the next three subsections ( H4.4II4.6) we provide three 
simulations generated by LUNA. The system configuration 
used in each case is discussed in the relevant subsections. In 
each case, we use a cadence of 1 minute and add on Gaus- 
sian noise of 250 ppm, in-lin e with the properties of K epler's 
short-cadence photometry (jKipping fc BakosH2011br ). Data 
are produced surrounding ±0.5 days of tb, for N consecu- 
tive transit epochs. 

These noised light curves are then fitted using a 
Metropolis-Hastings MCMC routine with 125,000 steps and 
a 20% burn-in to give 10 final points for building the pa- 
rameter posteriors. Jump-sizes are set to be equal to the 1- 
a uncertainties (determined through preliminary runs) and 
the final parameter values are given by the median of each 
parameter with uncertainties of 34.15% quantiles either side. 

For each case, we perform three fits with fixed assump- 
tions i) no moon is present i.e. s = and Ms/Mp — ii) 
the moon is prograde (i.e. isB is bounded to be in the range 
< isB ^ 7r) iii) the moon is retrograde (tt < isB ^ 2-7r). 
In reality, only one case is genuine and the two other act 
to show how distinguishable each scenario really is. We al- 
ways assume a circular orbit for the moon and planet for 
simplicity, which removes four parameters. 

In a l l cases , we utilize a powerful trick pointed out in 
iKippingl l|2010bT ). The detection of a planet-moon system al- 
lows one to determine the absolute mass and radius of the 
host star through Kepler's Third Law alone. Armed with 
these parameters, the planetary and lunar physical dimen- 
sions may also be obtained, thus replacing the traditional 
need for invoking spectroscopy and stellar evolution mod- 
els. The only additional information one requires is the ra- 
dial velocity semi-amplitude induced by the planet-moon 
barycentre, K. As the observational uncertainty on this pa- 
rameter scales with ~ 1/VArv, where N-rv is the number 
of RV measurements, it can be determined to high precision 
by simply repeatin g the rad i al veloc ity observations. 

As shown in Kipping (2010b), the absolute dimen- 
sions of the star have scalings M, ~ Kt and P„ ~ K t , 



and so [u(M,)f ~ [3a(K,)] 2 + [cr(AQ] 2 and KP*)] 2 ~ 
[o"(_ftT*)] 2 + [<r(P»)] 2 , where the dashed terms denote the 
value derived assuming no error on K t . In an example 
shown later (see ig3J, K* = 6m/s and a(M»)/M» ~ 15% 
and a(R' t ,)/R' t ~ 6%. This example is a singular case but 
nonetheless serves the purpose of illustrating a feasible situ- 
ation. In this example, the assumption that K* contributes 
negligible error into the stellar mass determination is valid 
for a(K»)/K» <2.5% (~ 0.15 m/s) and for the stellar ra- 
dius a(K,)/K, <7.0% (~ 0.4 m/s). These RV precisions are 
comparable to values being reported by those using high- 
preci s ion facilities such as Keck and HARPS l|Vogt et al.l 
|2010| ; lLovis et al .1 [20 1 lh . More generally, once an exomoon 
system is found, the uniqueness and importance of the sys- 
tem would certainly warrant the exploitation of such spec- 
trographic resources and so high levels of precision in K, 
can be expected. 

We therefore assume that the uncertainty on K* is much 
less than the uncertainty on the photometrically determined 
exomoon parameters in what follows. Th e inclusion of the 
determination of M* and R* using the iKippind (|2010bh 
method will demonstrate the feasibility of the technique. 

We note that s and Ms/Mp are positive definite 
and thus always yield a non-zero result. To identify cases 
where this bias i s crea ting a false-positive, we apply the 
iLucv fc Sweenevl i|l97ll ) test. If s or Ms/Mp yield a false- 
alarm-probability below 5%, then we quote the MCMC re- 
sult as usual. Otherwise, we quote the 95% confidence upper 
limit on those parameters. 



4.4 HZ-Neptune with a Close, Prograde Moon 
around an M2 Star 

4-4-1 Simulation and Fitting 

As a first example, we consider a Neptune-mass and radius 
planet in orbit an M2 star in the habitable- zone (K* = 
6.0ms -1 ). Using the stellar parameters from ICoxl |2000), 
we assign M» = 0.40 M© and P, = 0.50 P Q . We place the 
planet in an orbital period of Pb* = 46.0 days in a circular 
orbit with is* = 90°. Quadrati c limb darkening coefficients 
were generated using a [Kuruczl f2 006l style a tmosph ere fol- 
lowing the methodology of Kipping fc Bakosl (|2011bl ) , giving 
mi = 0.3542 and U2 = 0.3607. In total, we generate N = 6 
transit epochs spanning 0.85 years of continuous photome- 
try. 

For the exomoon, we consider an Earth-mass and radius 
moon on a close, prograde orbit around the host planet. We 
therefore select asB to be 5% of the Hill radius, correspond- 
ing to ~ 2 Pp. Through Kepler's Third Law, the orbital 
period is Psb = 0.3142 days. Such a moon would be ex- 
periencing large tidal dissipation and may not necessarily 
persist for Gyr, however, this is irrelevant for the purpose of 
generating some example transit light curves. We choose to 
place the moon in a slightly inclined orbit of isB = 92.93°, 
corresponding to bsB = —0.1. We also use fiss = 5° and 
assume a circular orbit. In the left panel of Figure [6] the 
simulated light curve is shown, before any noise is added. 

Due to the low impact parameter of the barycentre and 
the high coplanarity, TDV-TIP effects will be small and 
thus a determination of the sense of orbital motion will 
be highly challenging (see £12.21 for explanation) . However, 
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coplanarity is a reasonable c hoice as highly inclined moons 
have reduced Hill stability (|Donn ison 201C)j) and thus are 
less like ly a-priori. The same argument is true for eccentric 
moons (|Domingos et al.ll2006l ). 

4.4.2 Results 

As discussed in £14.31 we performed three fits for prograde, 
retrograde and no moon assumptions. A comparison of the 
fitted parameters for each assumption is presented in Ta- 
ble \E\ including the x 2 and Bayesian Informatio n Crite- 
rion (BIC) of the best- fits for each model. BIC |Schwar j 
Il978l : iLiddle et al.ll2007f ) is a tool for model selection which 
severely penalizes models for including more degrees of free- 
dom (i.e. Occam's razor) and is defined by: 

N / robs rmodcl \ 

BIC = ^(^ J +fclnA (44) 

where fi denotes flux, Oi is the associated uncertainty, 
N is the number of observations and k is the number of 
free parameter in the fit. We find that the prograde moon 
is only marginally preferable to the retrograde model, as 
expected. However, the no-moon model is clearly a poor fit 
and would be rejected. An F-test finds the prograde moon 
model accepted over no-moon model with a confidence of 
24.0-O", representing a very secure detection. 

The fit assuming a prograde orbit is shown in the right 
panel of Figure [6l where the thick black line represents the 
actual best-fit (i.e. not the original simulation light curve). 

We find that all of the parameters from the prograde fit 
are consistent with the true model. The physical parameters 
are generally poorly constrained, but apparently sufficient to 
identify the moon as rocky rather than icy. 



fit are consistent with the true model. The physical param- 
eters are generally well constrained, in contrast to the close, 
prograde moon considered in £14.41 This is due to osb/R* 
being significantly larger and thus determined to a higher 
precision which feeds into the other parameters. 

4.6 HZ-Neptune without a Moon around an M2 
Star 

4-6.1 Simulation and Fitting 

To complete the picture, we simulate the same case as the 
previous two subsections but removing the moon altogether. 
The data are then fitted assuming the three models as be- 
fore. The purpose of this data set is to show what happens 
when LUNA is implemented on control data and to ensure 
such data does not produce false positive moon detections. 

4.6.2 Results 

A comparison of the fitted parameters for each of our three 
model assumptions is presented in Table including the x 2 
and BIC values of the best-fits for each model. We find that 
the no-moon model gives the lowest BIC by a considerable 
margin, as expecteqj. The retrograde and prograde models 
produce non-sensical values for most parameters, especially 
the physical parameters, exploring various scenarios with 
no clear minimum in x 2 space. Note that s and Ms/Mp 
are positive definite and therefore are slightly skewed from 
a zero value, but not significantly so. An F-test finds the 
retrograde moon model (the best moon fit) accepted over 
no-moon model with a confidence of 2.0-ct, which indicates 
that the detection criterion is clearly in excess of 2-a. The 
use of BIC is therefore more accurate as a model selection 
tool for exomoon detection. 



4.5 HZ-Neptune with a Far, Retrograde Moon 
around an M2 Star 

4-5.1 Simulation and Fitting 

Our second example is identical to the previous one in £|4.4l 
except now we push the orbiting moon into a more dis- 
tant, retrograde orbit. We place the moon at 90% of the 
Hill radius (and thus i nside the 93.09% limit determined by 
iDomineos etail l|200rJ ll. 



4-5.2 Results 

A comparison of the fitted parameters for each of our three 
model assumptions is presented in Table [6l including the \ 2 
and BIC values of the best-fits for each model. We find that 
the retrograde moon is only marginally preferable to the 
prograde model, as expected. However, the no-moon model 
is clearly a poor fit and would be rejected. An F-test finds the 
retrograde moon model accepted over no-moon model with 
a confidence of 50.3-cr, representing a very secure detection. 

The fit assuming a retrograde orbit is shown in the right 
panel of Figure [3 where the thick black line represents the 
actual best-fit (i.e. not the original simulation light curve). 

We find that all of the parameters from the retrograde 



4.7 General Observations 

The simulations above are all for M2 stars with habitable- 
zone Neptune-like plane ts and Earth-like moon s. These pa- 
rameters were chosen as iKipping et all (2009c) have shown 
that such cases are optimally detectable for habitable-zone 
scenarios. We also tried the same configurations but using a 
K5 dwarf (M» = 0.67 Mq and R, = 0.72 R e ) but the longer 
period of the habitable-zone (half as many transits in the 
same time window) combined with lower radius and mass 
ratios for both the planet and moon meant we were unable 
to find convergent fits. This is not surprising an d echoes 
the motif of t he MEarth project (llrwin et al.ll2009T ) and the 
predictions in lKipping et alj (|2009cT ). 

In reviewing our fits, we find that for quantities such as 
Ms/Mp and s, which are positive definite, an overestimation 
of their value is common due to the boundary condition that 
they are greater than zero and the generally low signal-to- 
noise. This is similar t o the situation for orbi tal eccentricity 
in radial velocity fits (jLucv fc Sweenevlll97lh . 

3 Note that the \ 2 is actually worst for the no-moon model, but 
the BIC is lowest. This is because BIC heavily penalizes models 
for including more degrees of freedom, see Equation 1441 
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Figure 6. Left panel: Simulation from LUNA of a habitable- zone Neptune with a close, prograde, Earth-like moon for an M2 host star. 
Right panel: Noised data (circles) of 250ppm per minute overlaid with best fit from an MCMC routine (solid). 



Table 5. Comparison of parameter estimates from various model assumptions used in the fits. Data generated for a Neptune with a 
close moon around an M2 star. 



Parameter 


Truth 


Prograde 


Retrograde 


No Moon 


x 2 


8786.35 


8773.48 


8775.23 


9414.23 


BIC 




8891.31 


8893.07 


9468.62 


Fitted params. 



P B , [days] 
r B * [BJD - 2454000] 

P 2 [%] 

b B * 
f B , [s] 
P S b [days] 
4>SB [°] 
s 

osb/R* 

bsB 
USB [°] 

M s /Mp 



Physical params. 

M, [Mq] 0.40 

R* [Rq] 0.50 

M P [Mj] 0.054 

R P [Mj] 0.346 

M s [M e ] 1.00 

R S [R e ] 1.00 

PS [gem -3 ] 5.5 



46.000000 
956.00000 
0.5071 
0.00 
15882 
0.3142001 
40 
0.0183 
0.139 
-0.1 
5 

0.058 
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5082+ ' 0042 
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0.015 
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12 
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o ni+ 016 

u - ui -0.17 
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26 
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0.000 
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0.52 
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1.31 
1.03 
6.5 



0.66 
0.23 
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0.11 
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0.22 
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Figure 7. Left panel: Simulation from LUNA of a habitable- zone Neptune with a distant, retrograde, Earth-like moon for an M2 host 
star. Right panel: Noised data (circles) of 250 ppm per minute overlaid with best fit from an MCMC routine (solid). 



Table 6. Comparison of parameter estimates from various model assumptions used in the fits. Data generated for a Neptune with a 
distant moon around an M2 star. 



Parameter 


Truth 


Retrograde 


Prograde 


No Moon 


x 2 


8529.95 


8518.65 


8518.68 


11477.45 


BIC 




8636.48 


8636.51 


11531.83 


Fitted params. 



P B , [days] 
t b * [BJD - 2454000] 

P 2 [%} 

b B * 
f B , [s] 
P S b [days] 

4>sb [°] 

s 

a S B/R* 
bsB 

n SB [°] 
M s /Mp 



Physical params. 

M* [M Q ] 

R* [Rq] 
M P [Mj] 
Rp [Mj] 
M s [M e ] 
Rs [Re] 

PS [gem -3 ] 
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956.00000 
0.5071 
0.00 
15882 
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40 
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5 
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Figure 8. Left panel: Simulation from LUNA of a habitable- zone Neptune with no moon. Right panel: Noised data (circles) of 250ppm 
per minute overlaid with best fit from an MCMC routine (solid), from the starting point of assuming a retrograde moon is in orbit. 



Table 7. Comparison of parameter estimates from various model assumptions used in the fits. Data generated for a Neptune with no 
moon around an M2 star. 



Parameter 


Truth 


No Moon 


Prograde 


Retrograde 


x 2 


8936.23 


8930.96 


8928.42 


8921.57 


BIC 




8985.34 


9046.26 


9039.41 


Fitted params. 


P B , [days] 
t b „ [BJD - 2454000] 

P 2 [%] 
T B * [s] 


46.000000 
956.00000 
0.5071 
15882 


46.0000271°;— 
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■S080+ - 0035 
u.ouou_ 0026 
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U.OUDO_ 0041 
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s; 0.0178 
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4 +59 
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Physical params. 
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M P [Mj] 
R P [Mj] 

M s [M®] 
Rs [R@] 

PS [gem" 3 ] 


0.40 
0.50 
0.054 
0.346 
1.00 
1.00 
5.5 




4 Qn +162940 

5 4+32.0 
°- 4 -5.0 
f, 0+293.1 

Q rj— |- 22 . 2 
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19.41?^ 
l-4li 5 3 6 
lS-51^4 1 


, 91 +5876270 
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^•'-4.4 
4 (-+3043.3 

q q-(-82 . 2 
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7(! +82460 
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i 7+69.3 
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5 DISCUSSION & CONCLUSIONS 

5.1 Comparison to Previous Exomoon Light 
Curve Simulators 

Other works in the scientific literature have made use of al- 
gorithms to simulate exomoon transit light curves. To our 
knowledge, there e xists several previous u s es of such rou- 
tines. Figure 4 of ISartoretti fc Schneider! (|l999l '> presents 
four simulations of planet-moon combined transits, although 
details on the metho ds use to produce these simulations 
are sp arse. Similarly, ISzabo et al.l (|2006l ) and ISimon et al.l 
l|2009h present figures with simulated combined transits, al- 
though the focus is on a custom definition of transit timing in 
the former and the effects of the Rossiter-McLaughlin effect 
in the latter (i.e. few details on the me thods use to generate 
the li ght curves are provided). Finally. ISato fe A sada (2009, 
l2010h present perhaps the most detailed account of a method 
to generate planet-moon transit light curves but specifically 
negate the longitude of the ascending node, orbital inclina- 
tion of the barycentre, orbital eccentricity of the moon, limb 
darkening effects and use constant velocity approximations 
to model the planet/moon motion. 

In all cases, there is no reference to the iFewelll (|2006l ) 
solution and we therefore assume this was not utilized in 
those works. Therefore, all of the previous models must be 
at leas t parti ally numerical in nature since the solution of 
IFewelll (|2006l ) is the first instance of an analytic solution 
for the area of common overlap between three circles. In 
contrast, LUNA is fully analytic and we indeed find no ap- 
preciable compu tational time dif f erence between using LUNA 
versus the usual iMandel fc Ago! l|2002l ) routine. 

5.2 Distinguishing Mutual Events from Starspot 
Crossings 

Visual examination of the simulations presented in §4] in 
particular Figure [6] reveals that mutual events (i.e. when 
the planet and moon eclipse one another, during the stel- 
lar eclipse) bear a close resemblance to the morphology of 
a starspot crossing (i.e. when a planet eclipses a starspot 
on the surface of the star, during the stellar ecl ipse). Exam- 
ples o f star spot crossing even ts can be found in lRabus et alJ 
l|2009h and lPont et all (|2007h . This observation leads to the 
question as to how one could distinguish between a bona-fide 
moon and a starspot. 

Although a detailed study of this question remains be- 
yond the scope of this work, we suggest here that there are 
several tools at our disposal to make such a determination. 
Firstly, starspots are expected to co-rotate with the stel- 
lar rotation period, P ro t, which can be found by, for exam- 
ple, trac king the long-term fl ux variation of the target star 
(e.g. see lHenrv fc Winnl l2008). Application of a spot model 
to the light curve can test whether the spots are consis- 
tent with a rotation period determined through other means 
l|Silva-Valio fc LanzallioTlh . 

A second method w e can use is to employ the expres- 
sions of lKippinal ( |2010bh . which make use of the system dy- 
namics to measure the mass and radius of the host star for 
exomoon systems. If the determined mass and radius are 
inconsistent other methods, such as spectroscopy combined 
with stellar evolution modelling, then it is likely the puta- 
tive exomoon signal is in fact due to starspots. Whilst a 



dynamically-determined stellar mass and radius consistent 
with the spectroscopic value is not a proof of an exomoon, 
one may evaluate the false-alarm-probability of a starspot 
coincidentally inducing features which when modelled with 
LUNA cause precisely the correct stellar mass and radius to 
be determined. Also, a moon may induce transit timing ef- 
fects which give rise a derived lunar density which can be 
compared to expectation. 

A third method we suggest is that starspot crossing 
events will vary in amplitude when viewed chromatically, 
whereas the exomoon signal will have a chromatic variation 
much lower (due to atmospheric molecular absorption) and 
typically u ndetectable (unless u tilizing a telescope such as 
JWST, see iKipping et ai1l2009dj) . As an example, one could 
observe the target simultaneously in the visible and infrared 
wavelengths to test for this chromatic variation. 

5.3 Summary 

In this paper, we have presented a new algorithm called LUNA 
for modelling the transit light curves of a single planet with 
a single moon transiting a star. LUNA was designed from the 
outset to satisfy several key criteria: 

■ Analytic (absolutely no numerical components) 

■ Dynamic (inherently accounts for all timing effects) 

■ Limb darkening incorporated (including non-linear 
laws) 

■ All orbital elements accounted for (e.g. eccentricity, lon- 
gitude of the ascending node, etc) 

As a result of being both precise and analytic, LUNA is a 
highly potent weapon in exomoon detection. All of the pre- 
viously predicted observational consequences, such as TTV, 
TDV-V and TDV-TIP, are inherently built into the routine, 
plus previously unconsidered effects such as ingress/egress 
asymmetry. This is done by modelling the reflex motion of 
the planet due the moon at every instance, which is ul- 
timately responsible for all of the aforementioned timing 
effects. Transits and mutual eclipses are accounted for us- 
ing non-linear limb darkening laws meaning LUNA models all 
known observational consequences of exomoons for transit- 
ing systems. Further, physical parameters of the star, planet 
and moon may be derived utilizing the dynamical trick de- 
scribed in iKippind (|2010bl ). 

We have provided simulations of our new algorithm 
and re-fitted realistic noised data with MCMC methods to 
demonstrate the inverse retrieval of paramete r s. In o ne ex- 
ample, these fits demonstrate that the lKippTn e (2010b) tech- 
nique is capable of estimating the stellar mass and radius 
down to ~15% and ~5% respectively. For the exomoon this 
becomes ~12% and ~6% respectively. Whilst this is only 
one singular realization, the simulations lend credence to the 
prospect of characterizing the internal structure and compo- 
sition of Earth-mass bodies with current instrumentation. 
We stress that these examples are primarily to illustrate the 
effectiveness of our new algorithm rather than estimate the 
range of configurations which can be conceivably detected, 
which will be considered in future work. 

In this work, we have assumed the moon is small 
(Rs/R* <S 1) but future work will focus on relaxing this 
constraint so that gas giant binaries may be modelled as 
well. 
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APPENDIX A: PLANET MOON MOTION 
Al Nested Two-Body Model 

We will here derive the relative sky-projected motions of the 
planet and the moon, as used in the LUNA algorithm. Our 
goal is to find analytic expressions for the terms Sp,, Ss* 
and Sps, which were introduced in t|2.1l We will assume only 
a single moon in our model and a single planet. To predict 
the positions of the planet and moon at any instant in time, 
one must solve the three-body problem. Since no general, 
analytic solution exists for this problem, one must choose to 
work in a restricted case where analytic expressions can be 
employed. 

In this work, we will us e the nested two-body model 
presented in iKippinel (|2010bh . This considers the motion of 
the moon to be independent of the star in the reference frame 
of the planet-moon barycentre (dubbed the "inner frame"). 
Therefore, in the inner frame the motions are Keplerian. 
The barycentre is then considered to also maintain Keplerian 
motion around t he host s t ar in t he "outer frame" . 

As shown in iKippind <|2010bh , this model is an excellent 
approximation for clsb < 0.537?sr, where clsb is the moon's 
semi-major axis around the planet-moon barycentre and Rh 
is the Hill radius of the planet. Since all b ounded prograde 
satell ites must satisfy asB < 0.4895-Rs l|Domingos et alj 
2006), then the nested two-body model encompasses all pro- 
grade moons. We also point out that one may switch the 
module computing the motions fairly simply, but through- 
out this paper we will use the nested two-body model. 



A2 Inner Frame 

Let us consider the inner frame first, comprising of the planet 
and moon orbiting a common centre of mass. We define the 
moon to planet-moon barycentre distance, at any instant in 
time, as tsb- The subscript notation represents S to B or 
"satellite" to "barycentre" (where barycentre is understood 
to be the planet-moon barycentre). This notation scheme 
will be employed throughout this paper so that there exists 
no ambiguity as to what is being referred to. 



tsb = r S B = 



agg(l - e|g) 
1 + esB cos fsB 



(Al) 



Where a$B is the satellite to barycentre semi-major 
axis, fsB is the true anomaly of the satellite around the 
barycentre and ess is eccentricity of the satellite's orbit 
around the barycentre. 

One of the simplest reference frames in which to view 
the orbital motion is when the line connecting the apoapse 
and periapse is the f-axis and the orbit lies entirely within 
the x-y plane. Placing the barycentre at the focus lying at 
{asses, 0,0}, the Cartesian coordinates of the satellite are 
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Figure Al. Orbital elements of a two-body system. For a tran- 
siting planet-star system, i ~ 90° a nd the observer is located at 
{X, Y, Z} = {0, 0, +oo}. Figure from I Murray & Correld] \20lA) . 



defined by (lower-case symbols are used to denote this is 
the simple view ol the system not yet accounting for an 
observer's viewing angle): 



Armed with these definitions, the transformation from 
the simple frame of {x, y, z} to the observed frame of 
{X, Y, Z} may be written afl: 



Rss = Pz(nss)Px(tSfl)Pz(wss)rss 



(A5) 



After simplification, this yields the following final Carte- 
sian coordinates: 



Xsb = rsB [COS flss COs(u>SB + /ss) 

— sin iss shifts sin(u;sB + /sb)] 
Ysb = rsB [sin Q,sbcos(ujsb + /sb) 

+ sin isB cos Qsb sin(u;sB + /ss)] 
Zsb = rsB cos isB sin(ws S + /sb) (A6) 

Note that Rsb = |Rssj = jrssj = rss- The planet's 
reflex motion can be found by simply using: 



RpbMp + Rsb M s = 



(A7) 





( X S B\ 




rss = 


vsb 






\Z8B 





rss cos f sb 
rsB sin fsB 




(A2) 



This is highly analogous to the usual equations describ- 
ing a planet-star system. In the case of a planet-star, one 
starts with the analogous version of Equation IA2I and then 
performs three rotations to account for the viewing a ngle an 
observer has to the system (|Murrav fc Correiall2010l ). These 
three angles are the argument of periapsis, w, the orbital in- 
clination, i, and the longitude of the ascending node, Q,. The 
rotations are performed sequentially in a clockwise sense for 
the z-x-z axes respectively and are shown in Figure |A"T1 

The same scheme will be adopted here for the planet- 
moon inner frame. For the planet-star case, i is generally 
close to 7r/2 for transits and thus is a large rotation. How- 
ever, in our case, we choose to perturb the orbit from the 
coplanar condition. Therefore, we define a small angle for 
the inclination rotation, tss, and use isB = 7r/2 — iss 
so that isB = tt/2 corresponds to a coplanar case, in an 
analogous way to what is used for the planet-star case. 
The three anglefl we need to rotate by are cjss, lsb and 
n.<?B. lMurrav fc Correial l|201fj ) showed that these rotations 
may be concisely expres sed in matrix notation. Following 
iMurrav fc Corr eia (2010), we denote the following clockwise 
rotation matrices about the x and z axes respectively: 



A3 Outer Frame 

To obtain the overall motion, we now need to account for 
motion of the planet-moon barycentre around the star. Plac- 
ing the star at the origin, the barycentre-star separation, at 
any instant, is given by: 



as*(l — ( r; ' . > , . . 

t-b* = |r B , j = — — ^ y- (A8) 

1 + SB* COS /s* 

Where as* is the barycentre to star semi-major axis 
(usually just dubbed a), /s* is the true anomaly of the 
barycentre around the star (usually just dubbed /) and es* 
is eccentricity of the barycentre's orbit around the star (usu- 
ally just dubbed e). 

In exactly the same way as was done in the inner frame, 
the barycentre's position can be described in a simple frame 
which is then rotated to account for the viewing angle of the 
observer. The simple frame has the barycentre at Cartesian 
coordinates: 



(A9) 









( r B * cos /b*\ 






- 


rs, sin/s* 








V o ) 



This is then rotated in the same way as before: 



/l 

~Px{<t>)= coscj) — sin(j!> 
V sin 4> cos cj> 



(cos ( 
sin<; 




— sin (j> Q\ 
cos cf> 
1/ 



(A3) 



(A4) 



4 ujgg and flsB take the range — > 2n. We also allow tss to take 
the range — > 2n although only 7r radians are required to uniquely 
define every point in space. This is done to prevent boundary 
conditions when fitting data. 



Pz(SJfl,)Px(! B ,)Pz(^,)rB, 



(A10) 



The last rotation, about z by Qb* , has no bearing on the 
transit because the light curve is defined by the star-planet 
separation only, which is invariant about a £-axis rotation 
since the observer lies at z — +oo. For this reason, there is 
no need to include this final angle in practice: 



Px(is*)Pz(us,)rs, 



(All) 



5 Note that lower case Cartesian coordinates are used for inner 
frames and upper case for outer frames 
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To compute the Cartesian coordinates of the planet and 
moon (rather than the planet-moon barycentre), we need to 
combine the results of the inner frame derivation with the 
process of the outer frame. This is achieved by considering 
that: 



rs» 
rp, 



Rsb + r B * 



(A12) 
(A13) 



The conversion from rs, — > Rs* is performed by re- 
peating the rotations used to go from vb* — > Rs* (and 
similarly for rp,). This gives: 

Rs* = Pz(fi B ,)Px(iB,)PzK)(R Sfl + r B .) (A14) 
Rp, = Pz(fi B ,)Px(!p.)PzK,)(RpB + r Sa .) (A15) 

Evaluating, this gives the satellite's position as: 

X s * = r B * cos(/ B » + lj b *) 

+ r SB cos(/ss + ljsb) cos(^b» + Osb) 

— r S B sin i S B sin(/ SB + uj S b) sin(cj B » + SIsb) (A16) 
Ys* = r B * cosi B » sin(/s* + lj b *) 

— r S B sin i B * cos i S B sin( f S B + ojsb) 

+ r,sB cos i B * sin i,sB cos(cj b , + Qsb) sin(/ss + ljsb) 
+ r,sB cos is, sin(oj B , + Qsb) cos(/sb + ljsb) (A17) 
Z,s* = r s * sin is* sin(/ B , +lub*) 

+ r,sB cos i B * cos i,sB sin(/g B + ljsb) 

+ r,sB sin i B * sin isB cos(lu b * + Qsb) sin(/ss + ljsb) 

+ r S B sin i B * sin(o; B , + fiss) cos(/ss + ljsb) (A18) 

Evaluating for the planet, we find: 

Xp* = r B » cos(/ B * + ws*) 

— rps cos(/ss + cjsb) cos(ws* + fiss) 

+ r BB sin iss sin(/ss + w SB ) sin(o; B , + Qsb) (A19) 
Vp* = rs* cos is* sin(/s» + ujb*) 

+ r PB sin i B * cos 15s sin(/ SB + ljsb) 

— r PB cos is* sin isB cos(ljb* + &sb) sin(/g B + ljsb) 

— r PB cos i B * sin(o; B * + Osb) cos(/sb + wsb) (A20) 
Zp* = tb* sin is, sin(/s» + ujb*) 

— rps cos i B * cos isB s'm(fsB + ljsb) 

— rp B sin is* siniss cos(cjb» + Qsb) sin(/g B + ljsb) 

— r PB sin i B * sin(w B * + fiss) cos(/ss + ljsb) (A21) 



A4 Sky-Projected Distances 

The usual planetary transit light curve is completely de- 
scribed by the sky-projected distance between the planet 
and the star. In the same way, the light curve of a planet with 
a moon is described by the sky-projected distances, defined 
as S — \fX 2 + Y 2 /R* (units of the stellar radius are used). 
Using r' B * = r B */R*, r' SB = tsb/R* and r' PB = r PB /R*, 
the sky-projected moon-star distance is: 



S%* = \t'b* cos(/ B , + ujb*) + r'sB cos(cj b „ + Qsb) cos(/ SB + lj S b) 

, 1 2 

- r SB siniss sin(aj B , + SB )sin(/ SB + lj sb ) 

+ \r' B , cos i B * sin(/ B , + oj b * ) 

- r's B cosisB sin is, sin(/ SB + cj sb ) 

+ r's B cosis* siniss cos(w B , + Qsb) sin(/ss + ljsb) 

+ r's B cos is, sin(o;s, + fiss) cos(/ S s + ^ss)j 



(A22) 



And for the planet: 



S%* = \r' B * cos(/ s , + ujb*) - t'pb cos(cj s » + £l SB ) cos(/ ss + ljsb) 

1 2 

+ t'pb siniss sin(oJs, + &sb) sin(/ss + lj S b) 
+ \r' B * cos is, sin(/ s * + uj b * ) 

+ r'p B cos i S B sins, sin(/ss + ljsb) 

— r'p B cos is, siniss cos(^s* + Osb) sin(/ss + ljsb) 

1 2 

— r'pB cos is, sin(oJs, + &sb) cos(/ss + wss) 

(A23) 

Another important term is the separation between the 
planet and moon alone: 

c2 {X P * -Xs*) 2 + (Y P * -Ys*) 2 

^ M 



Ssp = [t'pb + r'sB? y [ cos(^s, + ^ss) cos(/ S s + lj S b) 



- siniss sin(cj s , + fiss) sin(/ S s + ^sb)J 
+ [ sin is, cos iss sin(/ss + ljsb) 

- cos is, siniss cos(cj s , + ^ss) sin(/ss + ljsb) 

- cosis, sin(o;s, + fiss) cos(/ S s + wss)] ' ^ (A24) 
A5 Repeating Units 

Inspection of the above equations reveals several repeating 
units. Substituting these units allows for a more simple and 
intuitive expression of the sky projected distances. The fol- 
lowing two units repeat in all three sky projected distances: 

P = r' SB [cos(^s, + S b) cos(/ S s + lj sb ) 

— siniss sin(o;s, + Oss) sin(/ss + w SB )| (A25) 

7 = r SB [cosis, siniss cos(cj b , + SB ) sin(/ S s + lj S b) 
+ cosis, sin(o;s, + fiss) cos(/ss + ljsb) 

— sin is, cos iss sin(/ss + wss)j (A26) 

Using these substitutions, the sky-projected distances 
are simplified to: 
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Ss* = [r' B * cos(/ B * + uj B *) + P] 2 

+ [r' Bt cos i s » sin(/ B * + ui B *) + "/] 2 (A27) 
S*p* = [r' Bt cos(/s* + ws.) - (M s /Ms)/3] 2 

+ [r' B , cosi B , sin(/ B * + w B *) - (M S /M P ) 7 ] 2 (A28) 

Ssp = [1 + (M s /Mp)] 2 [/3 2 + 7 2 ] (A29) 

The substitution terms may be split into terms which 
vary with time and those which do not, using a matrix I 
defined as the "inner" matrix: 



I = 



~p 




A 


/&" 




sin(/ ss + lo S b) 


.7 




71 


72 




C0s(/,SS + OJ.Sb)_ 



(A30) 



Where: 



/3i = — r' s sin iss sin(aj B * + Q,sb) (A31) 
Pi = r' s cos(cj b » + Q S b) (A32) 

71 = r' s [cos i s * sin i S s cos(o;s* + fiss) - sini B * cos i SB ] 

(A33) 

72 = cos is* sin(o;s* + 0,sb) (A34) 

As the matrix form provides a natural way to describe 
the motions, it is convenient to extend it to the barycentric 
motion as well: 



S%« = [O2 + Ii] 2 + [Oi + I 2 ] 2 



(A35) 



Sp* = [0 2 - (Ms/Mp)^] 2 + [d - (M s /Mp)I 2 ] 2 (A36) 
Where: 



01 = = r' Bt cos i B » sin(/s» + cj b ») + 

= tpi sin(/ B * + oj B *) + ip2 cos(/s» + w B „) (A37) 

02 = e = + r' B , cos(/s, + u) B *) 

= ei sin(/s* + oj b ,) + e 2 cos(/ B * + w B *) (A38) 
Where we have used: 



ipi = r B * cos is* 

</>2 = 
£1 = 
£2 = r B » 



(A39) 
(A40) 
(A41) 
(A42) 



Here, ipi = 6 B *, the impact parameter of the barycen- 
tre across the face of the star, and 62 = a Bt /R, for a cir- 
cular orbit. These two parameters are typically fitted for in 
planetary light curve analysis, but it is widely known that 
they share a strong correlation. In matrix form we have the 
"outer" matrix as: 



O 









0" 




sin(/ B * + oj B *) 









£2_ 




cos(/s» + CJB»)_ 



(A43) 



